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1. Introduction 

The original idea of the Thomson problem([TJ [2] was to find equilibrium positions of N 
charges constrained to a spherical surface interacting with Coulomb's law. From Earnshaw's 
theorem (3] |4] we know, however, that a collection of point charges cannot be maintained in 
a stable stationary equilibrium configuration solely by their mutual electrostatic interactions. 
Nonetheless, there is a steady interest to find minimum energy configurations of N charged 
particles on a sphere, see e.g. 16] 18] [TUl [TT] . There is also an interactive Java applet 
by Bowick et al. 02] from Syracuse university (NY) to find minimum configurations for the 
more general r~" potential using several different minimization algorithms. 

In this paper, we generalize the Thomson problem to arbitrary curved non-self- 
penetrating parametrized two-dimensional surfaces that are embedded in three-dimensional 
Cartesian space. We present the mathematical details to study the motion and minimum 
energy configurations of an arbitrary number of charged particles constrained to these 
surfaces, and we briefly describe how to implement the resulting N-body simulation using the 
Compute Unified Device Architecture (CUDA) of todays graphics boards. The generalization 
to r~" potentials is left as exercise for motivated students having basic experiences with 
electrostatics and differential geometry. 
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The structure of the paper is as follows. In Sec. [2] we briefly discuss the differential 
geometric details for the parametrization of trajectories on curved surfaces. The equations 
of motion for a free particle and of particles interacting by means of the Coulomb law are 
described in Sec. [3] In Sec. |4] we discuss some implementation details for the N-body 
simulation and the subsequent visualization. Several examples and some feasible exercises 
are presented in sections [5] and [6] 

The source code to reproduce the examples in this paper is written in C/C++/CUDA and 
is freely available from (http://www.vis.uni-stuttgcirt.de/~muelleta/ChaPaCS). It 
can be compiled on Linux and Windows systems. 



2. Particle trajectories on curved surfaces 

The trajectory 7(f) of a particle that is constrained to a curved surface is described in 
a straightforward manner by means of differential geometry. For that, consider a two- 
dimensional surface in the common three-dimensional Cartesian space M? parametrized by 
the function 

f:Rbl)^I 3 , (« 1 ,M 2 )i->f(M 1 ,M 2 )., (1) 

see figure [T] 




Figure 1. Parametrized surface f (U) Cll 



The intrinsic geometry of this surface is given by the first fundamental form 
df_ df_\ _/dt_ dt_\ _ 

which defines a metric (g) (/ = gy on f. Here, (•, •) denotes the standard scalar product in ]R 3 . 
The corresponding Christoffel symbols of the first kind 

2r _ d 8jk | dg ik dgjj 
' h du' duJ du k 
describe the effects of parallel transport in curved spacetime. Later on, we also need the 
Christoffel symbols of the second kind which follow from the Christoffel symbols of the first 
kind by raising the last index using the inverse metric = (g _1 ),-,■, 



^j=L g k " T ^n (4) 
m=l 

The inverse of the 2x2 matrix (g),-j = g[j reads 

1 f g22 —gu 
det(g) V -821 gu 



S-' = A^ l 22 ~ m ) (5) 
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with determinant det(g) = gugn- gngu- 

A particle trajectory on f can now be represented by the parametrization y : u(t) = 
(w 1 (t),u 2 (t)) £ U forf g [ti,tf], where t usually denotes time. The time derivative of f yields 
the tangent of y in R 3 , 



f= -f («(«)) = La^- 



dt 



du' dt 



The scalar product f 2 = (f,f) can be represented using the first fundamental form: 



• 2 / dt dt \ , i- 2 / dt dt \ 1 2 



df df 

(9m 2 ' (9m 2 



(" 2 ) 2 = E ^7"'"'' 



Several examples of surface parametrizations are listed in Appendix D 



(6) 



(7) 



3. Particle motion 



3.1. Equation of motion for free particles 

As long as there is only a single particle or there is no mutual interaction or external force, the 
trajectory of a single particle of mass m on the surface f follows from the Euler-Lagrangian 
(EL) equations 

d dh dh r , 

°-dtd¥-^ k = ^ (8) 

with Lagrangian L = yf 2 and f 2 from (|7ji. The derivative of the Lagrangian L with respect to 
u k and the subsequent time derivative yields the first part of the EL equations, 



d dh 
dt du k 



2 

E 



dr 



gik ) u' + ^gikii' 



For the partial derivative of L with respect to u k we obtain 
dh 



(9) 



(10) 



The derivatives of the metric gjj are detailed in Appendix B Equations (|9|l and ( 10 1 contracted 
with the inverse metric g'-i lead to the geodesic equation 



= u k + E rf ; -M'M ; ' (11) 
of a free particle constrained onto f. As expected, the mass m of the particle drops out. In 



flat space, g,j — 5 i; with 8,j being the Kronecker-5, the Christoffel symbols vanish and ( 1 1 
represents the parametrized form of the force-free Newtonian equation. 



3.2. Equation of motion for charged particles 

The Lagrangian for a charged particle q under the influence of the Coulomb interaction of 
other charged particles located at f„ reads 



N 



1 



L= m f 2 -Y 

2 ll f - f «l 



(12) 
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with the distance r .— ||f — f„| = \J (f(u) — f(u„) ,f(u) —f(u n )) between a particle at f(u) and 
a particle at f(u„). For the Euler-Lagrangian equations we need the partial derivative of the 
reciprocal distance 

—r^ = - <f-f n ,f-f„) = -r- 3 / f-f„,^\ . (13) 

duJ 2 duJ x ' \ dull 

As r~ 1 is independent of u J , the equation of motion of a charged particle under the influence 
of N other charged particles on f is given by 




3.3. Kinetic and field energy 

The total kinetic energy T of all charged particles is just the sum of the kinetic energies of the 
single particles: 

m N 

T = \Y$- (15) 

(=1 

The field energy W follows from 



3.4. External electric and magnetic fields 

The influence of an external electric or magnetic field onto a single charged particle on a 
curved surface is expressed by the Lagrangian 

L = |f 2 -#(f,0+<z(f,A(f,0) (17) 

with the electric potential 0, the magnetic vector potential A, and the relation E = — V(j> and 
B = V x A. The Euler-Lagrangian equations yield 

°-*+iw-<(*+ i **£y+<(£-wy- <i8) 

If the magnetic vector potential does not explicitly depend on time t, the last term of ( 18 1 can 
be dropped. 



4. N-body simulation and visualization 

As long as the number of particles in an N-body simulation is in the order of a few thousands, 
we do not need any specific acceleration algorithm and/or approximation procedure, but we 
can calculate the N-body interaction by brute force: every particle interacts with every other 
particle. The computation time, however, increases quadratically with the number of particles. 
A first step to accelerate the computation is to handle each particle by a separate compute unit 
that has to integrate the equation of motion ( fT4] i for this particle. The only prerequisite is 
that all compute units must have access to all particle positions and velocities which can 
be achieved using a shared memory system. Today, virtually all standard home computers 
and even high-end smartphones have at least a dual core processor inside that have access 
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to shared memory. Parallelization of computation can then be realized, for example, by the 
programming interface OpenMP\\3\ that splits the computation into several threads. 

Much higher parallelization can be achieved using the compute capability of modern 
graphics hardware. The Compute Unified Device Architecture (CUDA) or the Open Compute 
Language (OpenCL) offer a C-like programming interface in order to use the graphics 
processing units (GPUs) for general purpose computations. Even without sophisticated 
algorithms for an efficient memory access, GPU computation leads to an enormous speed- 
up. Together with the Open Graphics Library (OpenGL) we can explore physical simulations 
at interactive frame rates. 

The basic structure of our implementation is shown in figure [2] The basic block is 
the GLUT [ 14 1 main loop which acts on key strokes and mouse events, and which initiates 
rendering new frames. The "CUDA" block includes all surface descriptions (f, di,gij,T k j) 



scene defintions 
single_sphere.ini 




render 



Surface 



surface 



use FBO as 
surface texture 



Main program 



GLUT 
main loop 



Frame 
Buffer 
Object 



initiate 



time step 



render 



particles 



CUDA 
i 



surface descriptions 
— I cs_sphere.ini i 
— ' cs_torus.ini i 
— ' cs frustum.ini i 



calculate new positions 
and velocities 



f 

Vertex 
Buffer 
Object 



particle mapping 



Figure 2. Basic structure of the program. 



and calculates a single time step using either a standard Runge-Kutta second or fourth order 
method, see e.g. Press et al. |15|. To calculate a time step, it uses the particle positions u k 
and velocities u k stored within a Vertex Buffer Object (VBO) which can be directly accessed 
by CUDA and OpenGL. The new particle positions can then be rendered directly or they can 
be mapped onto their corresponding surfaces. This mapping is realized by means of a Frame 
Buffer Object (FBO). This FBO is an internal rectangular image (display) which, in our case, 
represents the domain U. This image is then used to texturize the surface. 

As an example, figure [3] shows particles (yellow splats) projected onto a rectangular 
image (FBO) that is used as texture for the sphere. The inner gray rectangle covers the whole 
domain of the sphere u 1 = (p £ [0, 2it), « 2 = i?£ (0, it). The dark red border slightly expands 
the domain to prevent particle splats near the domain's boundary from being clipped, see 
figure |4] The splats have to be distorted by means of the inverse metric of the surface to let 
them appear as circular splats when mapped onto the surface. 

The "Surface" block is responsible for drawing the surfaces themself. For that, we 
uniformly sample the domain using quads. These quads are then split into two triangles which 
are transformed by means of the surface function f of ([T]i within a so called vertex shader (see 
OpenGL Shading Language lfT6ll ). Figure [4] (right) shows the resulting wireframe of a sphere. 

To keep the code simple, we do not have any sophisticated scene description language 
but implement each scene in a separate ".inl"-file. A specific scene and its particular scene 
parameters have to be chosen at compile time. Each ".inl"-file must have three functions: 
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Figure 3. Particle rendering as Gaussian splats with radii d$ = ds/r and d(p = cfa/(rsin$) 
for fixed value ds. This FBO image can be mapped onto a sphere, see figure [4] The inner 
gray rectangle covers the whole domain of the sphere «' = <p 6 [0,27r), u 2 = i? 6 (0,tt). The 
dark red border (color online) prevents particle splats near the domain's boundary from being 
clipped. 




Figure 4. Particle rendering as Gaussian splats without (left) and with (center) extended 
domain, compare figure|3] Without extended domain the splats are clipped. Right: wireframe 
view of the sphere. 



init_Objects() defines all surfaces and assigns IDs to them; set_Particles () registers 
the number of particles, their parameters (mass, charge, initial position, initial velocity), and 
the ID of the surface they belong to; set ^Supplement () offers the possibility to change the 
camera parameters or the size of the window. The particle data could also be loaded from file 
when starting the program. 



5. Examples 

In the following examples, we use the explicit second-order Runge-Kutta method without 
step size control to integrate the equation of motion ( 14 1. The units in use are explained in 



Appendix A 



As a first example, we consider a Thomson problem situation where N = 128 charged 
particles (q = lOe) are located on a sphere of radius r = 1 mm. At the beginning of the 



Charged particles constrained to a curved surface 



7 



simulation, the particles are randomly distributed and have zero initial velocity, see figure [5] 
The field energy W\ t=0 /(Km e ) w 8.66684- 10 5 mm and the kinetic energy T = 0. To obtain 
a minimum energy configuration, we add a linear friction term with frictional constant tj = 10, 
see , 



Appendix C At simulation time (t = 4s), the field energy has dropped to W\t=4/(Km e ) 



7.39316 • 10 3 mm~ l and the kinetic energy has reduced to T/m e « 8.11308 ■ 10~ 3 mm 2 s~ 2 , 
which is nearly impossible to observe visually, however. After about 19.5 seconds, the kinetic 
energy has dropped below T/m e — 10~ l2 mm 2 s~ 2 . 



Figure 5. N = 128 charged particles on a sphere at simulation times t = {0, 0.2,1,4},; with 
step size At = 10~ 3 s and velocity-dependent friction 7] = 10 (see |Appendix C| . 

In table [T] we compare our minimum energy results for several number of particles N 
with the literature values. Glasser and Every's [17| estimation formula 

E (N) = y (l-aN- l/2 + bN- 3/2 ) (19) 

with parameters a = 1.10461 and b = 0.137 from ifTSll already gives a good approximation. 
The lowest energies for 1 10 < AT < 200 from Morris et al. [ 18 1 are determined using a genetic 
algorithm. 



N 



W 



lOOWiv 



100£(AO 



112 
128 
161 
200 



5.618044887 • 10 5 
7.393007455 • 10 5 
1.183308476- 10 6 
1.843885657 • 10 6 



5.618044882 • 10 5 
7.393007443 • 10 5 
1.183308474- 10 6 
1.843884272 • 10 6 



5.618079704 • 10 5 
7.392951914- 10 5 
1.183308683 • 10 6 
1.843881429 • 10 6 



Table 1. Minimum field energies W for N particles with charge q = lOe, time step At = 10~ 3 i' 
and frictional constant 77 = 10; compared with values Wn taken from Morris et al. 1181 and 
estimation given by j!9) . The factor 100 is because we use q = lOe instead of q = le. 

In the second example, we consider = 1024 particles (q = lOe) on a torus with 
radii R = 2 mm and r — 0.9 mm. As before, the particles are randomly distributed at the 
beginning of the simulation, see figure [6] At t ss 30s, the kinetic energy has dropped below 
T jm e = 10 _8 mm 2 s -2 and the field energy reads W /(Km e ) w 2.19652- 10 7 mmT 1 . Expectedly, 
all particles move to the outer side of the torus. Again, to stress the validity of our code, we 
compare our minimum energy configurations with the literature, see table [2] As can be seen, 
we are in good agreement with the literature values and sometimes we have found even lower 
energies. 

The dynamic evolution of the Af = 415 example of table [2] can be read from figure [7] 
At the beginning of the simulation, there is a very short peak of high kinetic energy due to 
the particles that are accelerated from rest. After less than 2 seconds, the kinetic energy has 
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Figure 6. N = 1024 charged particles on a torus at simulation times t = {0, 8}.v with step size 
At = 10~ 2 s and velocity-dependent friction 7] = 200. 



N 


a=R/r 


W 




IOOWjv 




20 


1.414 


1.029846718 


10 4 


1.029846689 


10 4 


20 


1.618 


1.098582566 


10 4 


1.098582529 


10 4 


100 


1.414 


3.082218338 


10 5 


3.082217005 


10 5 


100 


1.618 


3.295927689 


10 5 


3.296043624 


10 s 


415 


1.414 


5.660049622 


10 6 


5.660070457 


10 6 



Table 2. Minimum field energies W for N particles with charge q = 10e, time step At = 0.01 s 
and frictional constant 77 = 50 on a torus with aspect a = R/r and radius R = 1 ; compared with 
values Wn taken from 1121 . 



dropped below T jm e = 1 mm 2 s~ 2 . Then, at least from the visual impression, the particles 
do not move any more. However, there is still kinetic energy in the system which does not 
dissipate uniformly. Only after about 100 seconds, the particles slowly settle down and reach 
a minimum energy configuration (plateau in figure [7] right). However, depending on the 
initial random configuration, the simulation does not always reach the same minimum energy 
configuration. Then, the particles have to be either randomly distributed again or they have to 
be given a small jerk by pressing a key. 




20 40 60 80 100 120 140 160 180 
t (seconds) 



7 5.660 18e+06 - 
£ 5.66016e+06 - 
E 5.66014e+06 - 
& 5.660 12e+06 - 
^ 5.6601e+06 - 

OA 

c 5.66008e+06 - 

ig 5.66006e+06 - 

5.66004e+06 — 





20 40 60 80 100 120 140 160 180 

t (seconds) 



Figure 7. Kinetic energy (left) and field energy (right) depending on time for the torus 
simulation with N = 415 particles, aspect a = 1.414, step size A; = 0.01, and a frictional 
constant r/ = 50. The loss of kinetic energy is also due to the dissipative RK2 method. 
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A more intricate example is shown in figure [8] where two torii are intertwined and the 
particles have either the same or opposite charges. The corresponding field lines are shown 
in figure [9] They were started close to the charged particles in the direction of the outer 
surface normal and were integrated with a constant step size along the total electric field of 
all particles. Please note that the field lines end after n = 800 steps by default. And, as the 
object surface is no equipotential surface, field lines could also penetrate the object surface 
non-orthogonally. 




Figure 8. /V = 2 x 512 particles (q = ±10e) on two intertwined torii with radii R = 2 mm and 
r = 0.9 mm of either the same (left) or opposite (right) charge. 




6. Exercises 

In the following we give a few suggestions for possible exercises that could be done using the 
ChaPaCS source code. 

• Consider the Thomson problem situation where N = 128 charged particles (q = lOe) are 
located on a sphere. How does the minimum field energy vary with the size of the radius? 

• Adapt the set_Particles () method of the single sphere example such that one half of 
the particles have q = lOe and the other ones have only q = 5e. What happens? 
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• Determine the influence of the aspect in the torus example on the minimum field energy. 

• Find the parametrization of an ellipsoid and determine the metric coefficients gij and the 
Christoffel symbols T^. Expand the code to handle ellipsoids and find the minimum 
energy configurations. 

• Construct a new scene with a small sphere hovering above a plane. 

• Construct a new scene with four small spheres at the corners of a quad (quadrupole). 
Study the distribution of particles on the spheres and the overall field lines. 

• Determine a minimum energy configuration for = 128 charged particles (q = lOe) 
on the sphere and save it to disk. Restart the program with this configuration. Set the 
frictional constant to zero and add an external electric field. Slightly increase the electric 
field strength and discuss what happens. 



7. Summary 

In this work we have developed the equations of motion of N charged particles that are 
constrained to a curved two-dimensional surface but interact in three dimensions as usual. 
We have also given a short introduction how to implement the particle simulation using 
the compute capability of todays graphics boards. The source code (C/C++/CUDA) of the 
prototype implementation ChaPaCS is freely available and can be easily extended by other 
curved surfaces. 



Appendix A. Units 

For numerical computations we should know which order of magnitudes we have to deal with. 
If we use the electron mass m e and the electron charge e as basis units, we could rearrange 
(12) to 

L e 2 » QQ„ = M. 2 » QQ n 

m e 2 4 We „t 1 ||f-f n || 2 i-Jf-y ^ ' 

where M = m/m e and Q = q/e are dimensionless. If distances are given in millimeters and 
time is given in seconds, the constant K reads 

^=-^^^^0.253"^, (A.2) 
An£Qin e AKm e s L 

where m e « 9.109 ■ 10~ 31 kg « 5.686 ■ 10~ 18 eV ■ s 2 /mm 2 = 5.686 aeV ■ s 2 /mm 2 , e w 
1.602 • 10~ 19 C, Hq = 47T10- 7 Vs/(Am), and c = 299792458 m/s. Furthermore, we have 
Km e ~ 1.440 • 10~ 18 eV ■ mm — 1.440 aeV ■ mm. Then, the field energy ( 16 1, for example, 
reads 

Km, tik -Ml 
and has dimension one over length. 
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Appendix B. Derivatives of the metric 

The partial derivative of gy with respect to the coordinate u k is given by 

d d I df df \ _ / d 2 t df \ / df d 2 f 

dl? 8i} ~ d^Xdu^diJ / ~ \JuWu 1, J^/ + \^ , ~du J du~J 



(B.l) 



This can also be written using the Christoffel symbols of the first kind: 
d 

j u j i gij=^iKi + ^ik,i- (B.2) 
The time derivative of the first fundamental form reads 

~g ik = ^u" = (T M + T kn ,) u". (B.3) 
at du" v ' 

The partial derivative of the inverse metric g nm can be derived from d^gn = d u t (gingjmg"" 1 )- 
With the property gi„g tm = 8" 1 and some index juggling, we obtain 

■^8 =-{ r ik§ + r ,k8 )■ (B.4) 
Appendix C. Particle motion with friction 

Let the particle motion be damped by a velocity-dependent friction Fr = — h(v)^ with the 
frictional function h depending on the value of the velocity v = ||v|| and v = f . Then, the 
generalized friction R% reads 

. , v dr , , . v d\ h(v) I ■ df \ ,„ , s 

R k = -h(v)- ■ ^- = -/i v - ■ ^- = — ^ (f, j- k ) (C.l) 
v ou k v ou k v \ du K I 

For the linear friction h(v) = T}v with frictional constant T] and the time derivative of f, 
compare d5), we obtain 

2 / df df \ 2 



The Euler-Lagrange equation (|8]l now reads 
d dL dL 

dt du k du k 



d dL dL 

0=-n^7v-^- R k, (C3) 



and the equation of motion ( 1 1 1 is given by 



2 

o = a k + 77 a k + £ r'ljiiui . (c.4) 



A friction that is quadratic in the velocity yields the term TjM*||f|| with ||f|| = yYdj=\ giju'ti-i. 

The unit of 77 depends on the unit of the coordinates u k . As thes e could be diffe rent like in 
the frustum case, where [u l = (j)] = rad and [u 2 — z] — length (see Appendix D.4 1, we should 
use two separate frictional constants. But here, we use the same numerical value and use the 
respective unit if necessary. 
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Appendix D. Surface examples 

The following surface parametrizations are given in standard form, which means that the 
components of f are with respect to the global Cartesian coordinate system of M?. Thus, 
f = / 1 ei+/ 2 e 2 + / 3 e 3 with 

ei = (l,0,0) r , e 2 = (0,l,0) r , e 3 = (0,0,l) r , (D.l) 

and the center of the object equals the origin of the global coordinate system. However, for 
more elaborate scenes, the basis {ei,e2,e3} can be oriented and translated arbitrarily. 



Appendix D.l. Plane 

A plane is definitely the most simple surface and can be parametrized by the coordinates 
(u 1 = x,u 2 — y) which, in general, are bound to a fixed domain x G [x ai xi,\ and y e [ja,^]- 
The surface reads f(x,y) — (x,y,0) T , and the metric is given by gn = gri = 1 and gn = 0. 
The Christoffel symbols of the second kind vanish identically. 



Appendix D. 2. Sphere 

The surface of a sphere with radius r can be parametrized by spherical coordinates (u 1 = 

<p,u 2 = &): 

(sin#cos<p \ 
sin # sin <p (D.2) 
COS# J 

with <p e [0,27c), # G (0,7c), and partial derivatives 

/ cos # cos <p \ j)f [ — sin # ship 
-r — =r\ cos # sine) , -r— = r sintfcosffl 
V -sin/ ) d( P { 

The metric coefficients read 

gn=r 2 sin 2 #, g n = 0, g 22 = r 2 , (D.3) 
and the only non-vanishing Christoffel symbols of the second kind are 

r| 2 = cot#, r 2 ! = -sintfcostf. (D.4) 



Appendix D.3. Ellipsoid 

While the parametric representation of an ellipsoid is similar to the representation of a 
sphere, the metric coefficients are quite cumbersome. Again, we use spherical coordinates 
(u 1 = <p, u 2 = i?) but now -& denotes the latitude angle, 

(acos#cosa> \ 
£cos#sin<p . (D.5) 
csintf J 

The partial derivatives are straightforward: 

/ — asin#cos<p \ fif { — a cost? sin (j? 
- — = — frsin#sinc> , = /?cos#cose> 

I ccostf d( P \ 
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The metric coefficients are straightforward and the Christoffel symbols of the second kind 
read 

j (b 2 -a 2 )c 2 cos(<p)sin(<p)cos(#) 2 2 a 2 b 2 cos(#) sin(#) 
Fn = p ' Tn= p ' 

. _ (fr 2 ~a 2 )c 2 cos((p)sin((p) , _ 

1 22 — , i 12 — — C0t», 

2 _ [(£ 2 -a 2 )c 2 sin(<p) 2 -£ 2 c 2 W£ 2 ]cos(#)sin(#) 

where p = a 2 fr 2 sin(#) 2 + [(a 2 -b 2 )c 2 sin(<p) 2 + b 2 c 2 ] cos(#) 2 . 
Appendix DA. Frustum 

A frustum is defined by the upper and lower radii (n,r 2 ) and the height h. With parameters 
u 1 = <p g [0,2jt) and u 2 =zE [0,h], we have 

/p(z)coscp\ 

f((p,z) = p(z)sin<p with p(z)=n — z. (D.6) 

The partial derivatives are straightforward: 

df I -p(z)sin(p\ ^ f / p'(z) cos <p 
dj=[ P(Z)C ° S(P J' ^=^'« si ^ 

where p'(z) = dp/dz = —{r\ — r 2 )/h. From these partial derivatives we obtain the surface 
normal 



1 



cos <p 



n ^' z ) = ^^^f^ sin( p • (D - 7) 

Vl+P'W 2 V -P'(z) 

The metric coefficients read 

gn = p(z) 2 , §12 = 0, g 22 =p'(z) 2 + l, 
and the non-vanishing Christoffel symbols of the second kind are 

{r\-r 2 ) 2 + hri(r 2 -r i ) , r 2 -r x 



Ml- 7T , u2 ' 1 12 



(r l -r 2 ) 2 +h 2 ' 12 (r2-n)z + An' 



Appendix D. 5. Torus 



A torus is defined by two radii, where R is the radius of the main circle and r is the 'thickness' 
radius, 

((/? + rcos#)cos<p \ 
(fl + rcos#)sin<p . (D.8) 
rsin# J 

Here, m 1 = # S [0,2n) and u 2 = (p g [0,27c), and the partial derivatives read 

/ — sin#cos<p \ / — (R + rcos&)sin(p 

^— = r I — sin#sin<p J'^ = ( (/? + >"cos#)cos(j!> 
f 1 \ cos & / \ 
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The metric coefficients are straightforward 

gn=r 2 , £i2 = 0, g 2 2 = (R + rcos$) 2 , (D.9) 
and the non-vanishing Christoffel symbols of the second kind read 
2 _ rsinfl . _ (R + r cost?) sin g 

1 12 — — IT"! q i 1 22 — • (1J.1U) 

/? + rcos?> r 
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